Generating facies probablity cubes

ABSTRACT

A method for generating one or more geological models for oil field exploration. The method includes receiving one or more well facies logs, a vertical facies proportion curve, a lateral proportion map, a variogram model and a global target histogram. The method then includes generating a facies probability cube using a modified Sequential Gaussian Simulation (SGSIM) algorithm, the well facies logs, the vertical facies proportion curve, the lateral proportion map and the variogram model. After generating the facies probability cube, the method includes matching the facies probability cube to the global histogram and generating the geological models based on the matched facies probability cube.

RELATED APPLICATIONS

This application is a continuation of and claims priority to U.S. patentapplication Ser. No. 12/837,936 filed 16 Jul. 2010, titled GENERATINGFACIES PROBABILITY CUBES, which in turn claims priority to U.S.provisional patent application Ser. No. 61/315,098 filed Mar. 18, 2010,titled METHOD FOR RESERVOIR MODELING USING A GEOSTATISTICAL APPROACH TOGENERATE 3D PROBABILITY CUBES OF BINARY FACIES WITH GEOLOGICALCONSTRAINT, both of which are incorporated herein by reference.

BACKGROUND

1. Field of the Invention

Implementations of various technologies described herein generallyrelate to techniques for modeling geological properties of the earthand, more particularly, to techniques for generating facies probabilitycubes that can be used with multipoint statistics to create a reservoirmodel.

2. Description of the Related Art

The following descriptions and examples are not admitted to be prior artby virtue of their inclusion within this section.

Geological modeling and reservoir characterization provide quantitative3D reservoir models based on available reservoir measurements, such aswell log interpretations, experimental results from core analysis,seismic survey and dynamic fluid flow responses from field observations(e.g., historic production data) and pressure change data. One type ofreservoir characterization or modeling technique is stochastic reservoirmodeling.

Stochastic reservoir modeling has gained popularity in modern reservoirmodeling software because of its ability to constrain its model based ona variety of reservoir data and its computational efficiency ingenerating complex reservoir models with millions of voxels. Stochasticreservoir modeling also allows users to quantitatively evaluateuncertainties in the model due to the lack of knowledge of thereservoirs. The data used to constrain the reservoir models instochastic reservoir modeling are primarily classified into twocategories: “hard data” or “soft data.” Hard data includes data such asthose measured in wells (e.g., well log data), which are considered tobe accurate information and should be honored during simulations. Softdata are not as accurate as hard data but typically have larger orbetter coverage of the reservoir. Facies probability cubes areconsidered to be soft data that have been derived from seismicattributes using well to seismic calibrations. These types of soft dataare important in guiding inter-well facies prediction and thus, they maybe used to reduce uncertainties in the decision making process forreservoir management.

Geostatistics provide variety of algorithms and tools to buildstochastic reservoir models. Generally, there are two approaches forusing geostatistics to build stochastic reservoir models: a pixel-basedapproach and an object-based approach. The pixel-based approach proceedsby gridding the reservoir into pixels (voxels) and simulating each pixel(voxel) in a random sequence. Unlike the pixel-based approach, theobject-based approach directly drops the facies (geobody) objects intothe simulated reservoir according to the specified geometric informationof these geobodies. The pixel-based approach provides increasedflexibility for constraining the model according to reservoir data butit often has poor shape reproduction in the final reservoir models. Incontrast, the object-based approach tends to generate more realisticshapes of geobodies; however, it becomes more difficult to constrain themodels according to local reservoir observations, particularly whenthere are dense well locations.

In the pixel-based approach, a sequential indicator simulation (SIS) isoften used to create facies models. However, a newly emergingpixel-based approach named Multi-point statistics (MPS) is gaining moreattention from modelers and is considered to be part of an advancedfacies modeling suite. MPS uses 1D, 2D or 3D “training images” asquantitative templates to model subsurface property fields. MPS modelingcaptures geological structures from training images and anchors them todata locations. As such, MPS takes advantage of a 2-point(variogram-based) geostatistical approach and an object-based approachto create flexibility in data conditioning while producing morerealistic shapes from the training images. MPS can then integrate softdata, such as a facies probability cube, to generate geological orreservoir models. The resulting geological or reservoir models can thenbe used for oil field explorations by identifying hydrocarbon depositsin the Earth.

In addition to being used in multipoint statistics facies modeling,probability cubes are also used for other modeling approaches, such asobject-based modeling and pixel-based Sequential Indicator Simulation(SIS), to further constrain the simulated earth models. As such, faciesprobability cubes play an important role in geological modeling orreservoir characterization. In particular, facies probability cubes canassist in geological modeling or reservoir characterization when welldata is scarce. However, automatically generating facies probabilitycubes that are geologically sound remains a challenge.

SUMMARY

Described herein are implementations of various technologies forgenerating facies probability cubes and using the facies probabilitycubes to generate geological models for oil field exploration. In oneimplementation, a method generating geological models may includereceiving one or more well facies logs, a vertical facies proportioncurve, a lateral proportion map, a variogram model and a global targethistogram. The method may then include generating a facies probabilitycube using a modified Sequential Gaussian Simulation (SGSIM) algorithm,the well facies logs, the vertical facies proportion curve, the lateralproportion map and the variogram model. After generating the faciesprobability cube, the method may match the facies probability cube tothe global histogram and generate the geological models based on thematched facies probability cube.

The claimed subject matter is not limited to implementations that solveany or all of the noted disadvantages. Further, the summary section isprovided to introduce a selection of concepts in a simplified form thatare further described below in the detailed description section. Thesummary section is not intended to identify key or essential features ofthe claimed subject matter, nor is it intended to be used to limit thescope of the claimed subject matter.

BRIEF DESCRIPTION OF THE DRAWINGS

Implementations of various technologies will hereafter be described withreference to the accompanying drawings. It should be understood,however, that the accompanying drawings illustrate only the variousimplementations described herein and are not meant to limit the scope ofvarious technologies described herein.

FIG. 1 illustrates a schematic diagram of a logging apparatus inaccordance with implementations of various techniques described herein.

FIG. 2 illustrates a flow diagram of a method for generating a faciesprobability cube in accordance with implementations of varioustechniques described herein.

FIG. 3 illustrates a flow diagram of a method for generating a multiplefacies model constrained by a facies probability cube in accordance withimplementations of various techniques described herein.

FIG. 4 illustrates a flow diagram of a method for generating a multiplefacies model constrained by a facies probability cube in accordance withimplementations of various techniques described herein.

FIG. 5 illustrates a flow diagram of a method for performing a modifiedSequential Gaussian Simulation (SGSIM) algorithm in accordance withimplementations of various techniques described herein.

FIG. 6 illustrates a computer network into which implementations ofvarious technologies described herein may be implemented.

FIG. 7 illustrates an example of a vertical proportion curve inaccordance with implementations of various techniques described herein.

FIG. 8 illustrates an example of a lateral proportion map in accordancewith implementations of various techniques described herein.

FIG. 9 illustrates an example of well facies data in accordance withimplementations of various techniques described herein.

FIG. 10 illustrates an example of a fluvial sand probability cube inaccordance with implementations of various techniques described herein.

FIG. 11 illustrates an example of two Probability Density Function (PDF)curves of a Beta distribution for different parameter controlledcontrasts in accordance with implementations of various techniquesdescribed herein.

FIG. 12 illustrates an example of a Cumulative Density Function (CDF)curve in an original probability space and a CDF curve in a normal scorespace in accordance with implementations of various techniques describedherein.

FIG. 13 illustrates an example of four facies training image inaccordance with implementations of various techniques described herein.

FIG. 14 illustrates an example of a reservoir model with four facies inaccordance with implementations of various techniques described herein.

FIG. 15 illustrates a diagram for creating a reservoir model inaccordance with implementations of various techniques described herein.

DETAILED DESCRIPTION

The discussion below is directed to certain specific implementations. Itis to be understood that the discussion below is only for the purpose ofenabling a person with ordinary skill in the art to make and use anysubject matter defined now or later by the patent “claims” found in anyissued patent herein.

The following provides a brief description of various technologies andtechniques for generating facies probability cubes. In oneimplementation, a computer application may receive a vertical faciesproportion curve, a lateral proportion map, a variogram model and aglobal target histogram to generate the facies probability cube. Thevertical proportion curve may indicate the amount in which the faciesdeposits change within varying depths. The lateral proportion map may bea 2 dimensional map in the XY plane that represents the variation oflateral proportions of certain facies in terms of an average proportionalong the Z direction. The variogram model may be a quantitative tool inconventional 2-point geostatistics that measures the spatial variabilityof a geo-variable. Generally, the vertical facies proportion curve andthe lateral proportion may be determined based on well facies data inwell facies logs. The variogram model may be determined based on anautomatic or computerized interpretation of the well facies data or itmay be received from a user based on his/her interpretation of the wellfacies data or other relevant data. The global target histogram may bereceived from a user and may be defined as probability distribution suchas a U-shaped beta distribution.

After receiving the vertical facies proportion curve, the lateralproportion map, the variogram model and the global target histogram, thecomputer application may generate a facies probability cube using amodified Sequential Gaussian Simulation (SGSIM) algorithm. The modifiedSGSIM algorithm may use the vertical facies proportion curve, thelateral proportion map and the variogram model to estimate faciesproportions in a facies probability cube. The modified SGSIM algorithmmay estimate the facies probability cube by increasing or reducing akriging mean based on the vertical facies proportion curve (VPC) and thelateral proportion map (LPM). The adjustment with regard to the verticalfacies proportion curve (VPC) and the lateral proportion map (LPM) maybe performed based on the following formula:

${m_{sk}^{*{new}}\left( {i,j,k} \right)} = {{m_{sk}^{*}\left( {i,j,k} \right)} + {\frac{\lambda}{1 - \lambda}\left( \frac{{V\; P\; {C(k)}} - {p^{*}(k)}}{\sigma_{T}} \right)} + {\frac{\lambda}{1 - \lambda}\left( \frac{{L\; P\; {M\left( {i,j} \right)}} - {p^{*}\left( {i,j} \right)}}{\sigma_{T}} \right)}}$

where m_(sk)*(i, j, k) is kriging mean at pixel (i, j, k), p(k)* isrunning proportion of the facies at the k-th Z-layer, VPC(k) isproportion read at the k-th Z-layer from the VPC, λ is a servo factor in[0, 1] that indicates the correction strength, m_(sk)*^(new)(i, j, k) isadjusted kriging mean, σ_(T) is the standard deviation of the globaltarget histogram, i.e., σ_(T)=√{square root over (c×p×(1−p))} andLPM(i,j) is lateral proportion at pixel location (i, j), which can beread from the LPM, and p*(i,j) is a running facies proportion calculatedfrom the column at (i, j). After generating the facies probability cube,the computer application may perform an inverse normal score transformon the facies probability cube in a Gaussian space to match the globaltarget histogram.

FIGS. 1-15 illustrate one or more implementations of various techniquesdescribed herein in more detail.

FIG. 1 illustrates a schematic diagram of a logging apparatus inaccordance with implementations of various techniques described herein.FIG. 1 shows a borehole 32 that has been drilled in formations 31 withdrilling equipment, and typically, using drilling fluid or mud thatresults in a mudcake represented at 35. A logging device 100 is shown,and can be used in connection with various implementations describedherein. The logging device 100 may be suspended in the borehole 32 on anarmored multiconductor cable 33. Known depth gauge apparatus (not shown)is provided to measure cable displacement over a sheave wheel (notshown) and thus the depth of logging device 100 in the borehole 32.Circuitry 51, represents control and communication circuitry for theinvestigating apparatus. Although circuitry 51 is shown at the surface,portions thereof may typically be downhole. Also shown at the surfaceare processor 50 and recorder 90. Although logging device 100illustrated herein is shown to be a wireline logging tool, it should benoted that other tools such as a logging while drilling tool may be usedin connection with various implementations described herein.

The logging device 100 may represent any type of logging device thattakes measurements from which formation characteristics can bedetermined, for example, by solving complex inverse problems. Thelogging device 100 may be an electrical type of logging device(including devices such as resistivity, induction, and electromagneticpropagation devices), a nuclear logging device, a sonic logging device,or a fluid sampling logging device, or combinations thereof. Variousdevices may be combined in a tool string and/or used during separatelogging runs. Also, measurements may be taken during drilling and/ortripping and/or sliding. Examples of the types of formationcharacteristics that can be determined using these types of devicesinclude: determination, from deep three-dimensional electromagneticmeasurements, of distance and direction to faults or deposits, such assalt domes or hydrocarbons; determination, from acoustic shear and/orcompressional wave speeds and/or wave attenuations, of formationporosity, permeability, and/or lithology; determination of formationanisotropy from electromagnetic and/or acoustic measurements;determination, from attenuation and frequency of a rod or platevibrating in a fluid, of formation fluid viscosity and/or density;determination, from resistivity and/or nuclear magnetic resonance (NMR)measurements, of formation water saturation and/or permeability;determination, from count rates of gamma rays and/or neutrons at spaceddetectors, of formation porosity and/or density; and determination, fromelectromagnetic, acoustic and/or nuclear measurements, of formation bedthickness.

In one implementation, the measurements obtained by logging device 100may include well facies data in the well logs (facies logs). Facies logsmay indicate the absolute presence or absence of a targeted facies atdifferent spatial locations. Targeted facies may include channels,levees, crevasses, shale background and the like.

FIG. 2 illustrates a flow diagram of a method 200 for generating afacies probability cube in accordance with implementations of varioustechniques described herein. In one implementation, method 200 may beperformed by a computer application. The following description of method200 is made with reference to method 500 of FIG. 5 and diagram 1500 ofFIG. 15. It should be understood that while method 200 indicates aparticular order of execution of the operations, in someimplementations, certain portions of the operations might be executed ina different order.

At step 210, the computer application may receive well facies data fromwell facies logs. Well facies logs may be well logs that have beenacquired at various locations in a survey area such as boreholes and thelike, as described above with reference to FIG. 1. Well facies logs mayindicate the absolute presence or absence of a targeted facies atdifferent spatial locations. An example of well facies data is providedin FIG. 9.

At step 220, the computer application may receive a vertical faciesproportion curve (VPC). FIG. 7 illustrates an example of a VPC. FIG. 15illustrates how the VPC 1510 may be used in method 200. The VPC mayillustrate the change of facies deposits in the earth with the variationof depths. The change of facies deposits in the earth may occur due to ageological sedimentation process, which may evolve periodically asprogradation and retrogradation patterns. These patterns may result whensea levels rise and drop alternatively. This phenomenon is ubiquitous insequence stratigraphy at large scales, but the cyclical nature of thepatterns may also be found at small scales. As such, the amount ofspecific facies deposits changes with the variation of geological timesor depths leading to a systematic vertical trend of facies proportions.For example, in a fluvial-dominated deltaic environment, the sanddeposit may be high at the bottom of a reservoir unit, become lower inthe middle of the reservoir unit and becomes high again at the top ofthe reservoir unit. Coarsening or fining upwards characteristics ingeological deposition for a certain facies is also a trend indicator ofthe deposition of that facies.

In one implementation, VPC can be calculated from well facies logs (fromstep 210) interpretation, geological conceptual models or analogssimilar to the reservoir under study. If well facies data are available,the well facies data may be combined with geological interpretations ofthe seismic area to generated interpreted well facies data. Theinterpreted well facies data may then be used as anchors in generatingfacies probability cubes, i.e. they indicate the absolutepresences/absences of the targeted facies at different spatiallocations: probability is 1 for its presence and 0 for its absence. Thevertical facies proportion curve can be calculated by extracting thefacies proportion along each vertical layer (i.e., Z-layer) of amodeling grid. In case of very sparse wells, the calculated proportioncurve may not be reliable and hence, its modification and editing may beperformed based on geological interpretations or analogs similar to thestudied reservoir.

At step 230, the computer application may receive a lateral faciesproportion map (LPM). FIG. 8 illustrates an example of a LPM. FIG. 15illustrates how the LPM 1520 may be used in method 200. The lateralproportion map may be a 2 dimensional map in the XY plane thatrepresents the variation of lateral proportions of certain facies interms of an average proportion along the Z direction. In this manner,the lateral proportion map may illustrate lateral trends such asprogradation or transition. In one implementation, the computerapplication may obtain the lateral facies proportion map byinterpolating facies proportion from each well facies log or bycalibrating seismic interpretations of each well facies log. Theinformation included in the lateral facies proportion map may be usefulwhen facies depositions have apparent lateral trends, such asprogradation or a transition as observed from received seismic data,well facies log interpretations or conceptual models.

The computer application may create lateral proportion maps based onwell facies data using any smooth interpolator, such as a krigingtechnique in 2-point geostatistics or moving averaging algorithms. Forinstance, at each well location, the computer application may calculatea facies proportion. The facies proportion value may then be consideredto be one of the hard data used as an anchor point to control the smoothinterpolator. In one implementation, the computer application may alsoadd interpreted trends into the LPM by adding data from “pseudo” wells.In any case, by integrating LPM into reservoir modeling, the computerapplication may be able to generate a facies probability cube that hasmore realistic results.

At step 240, the computer application may receive a variogram model froma user. The variogram model may be a quantitative tool in conventional2-point geostatistics that measures the spatial variability of ageo-variable. The variogram model may be used to control the anisotropyof the spatial distribution of a geo-variable by changing thecorrelation ranges and the orientations of major/minor axis of thevariogram model. As such, the variogram model may be specified by a userto govern the anisotropic distribution of probability values in thefacies probability cube.

The variogram model may be derived by the computer application based onwell facies data or may be specified by a user based on his/herinterpretation of the well facies data or other relevant data. If thevariogram model is derived by the computer application, the computerapplication may derive the variogram model from the vertical faciesproportion curve and the lateral facies proportion map. The verticalfacies proportion curve may be used to infer the vertical range of thevariogram model, and the lateral facies proportion map may be used toinfer the horizontal correlations of the variogram model. If a localazimuth field that reflects the orientation of facies deposition isavailable, the computer application may also use the local azimuth fieldto build the variogram model.

Generally, the computer application may receive one of four commonlyused variogram models: spherical, exponential, Gaussian and powervariograms. It has been shown that the Gaussian variogram model, ascompared with the other three variogram models, may create smootherresults in generating the facies probability cube and therefore may bemore suitable to generate facies probability cubes. However, it shouldbe noted that method 200 described herein is not limited to the Gaussianvariogram model.

At step 250, the computer application may receive a global targethistogram from a user. The global target histogram may be a probabilitydistribution used to constrain the facies probability cube. In oneimplementation, the global target histogram may be a U-shaped betadistribution. As such, the two end peaks of the U-shaped betadistribution may correspond to the likelihood of facies presence/absencein the entire probability cube. Additionally, the mean of the Betadistribution may be forced to match the global facies proportion in thereservoir while the variance of the Beta distribution may be auser-defined parameter that allows the user to control the variation or“sharpness” of the resulting probability cube.

The probability density function (PDF) of Beta distribution is definedin classic statistics as:

$\begin{matrix}{{f\left( {x,\alpha,\beta} \right)} = \frac{{x^{\alpha - 1}\left( {1 - x} \right)}^{\beta - 1}}{\int_{0}^{1}{{u^{\alpha - 1}\left( {1 - u} \right)}^{\beta - 1}\ {u}}}} \\{= {\frac{\Gamma \left( {\alpha + \beta} \right)}{{\Gamma (\alpha)}{\Gamma (\beta)}}{x^{\alpha - 1}\left( {1 - x} \right)}^{\beta - 1}}} \\{= {\frac{1}{B\left( {\alpha,\beta} \right)}{x^{\alpha - 1}\left( {1 - x} \right)}^{\beta - 1}}}\end{matrix}$

where Γ is the gamma function. The expected value, second moment andvariance of a Beta random variable X with parameters α and β are:

${E(X)} = \frac{\alpha}{\alpha + \beta}$${E\left( X^{2} \right)} = \frac{\alpha \left( {\alpha + 1} \right)}{\left( {\alpha + \beta} \right)\left( {\alpha + \beta + 1} \right)}$${{Var}(X)} = \frac{\alpha \; \beta}{\left( {\alpha + \beta} \right)^{2}\left( {\alpha + \beta + 1} \right)}$

To match a specified global facies proportion (p) and a user-definedtuning factor (c), the following relationships may be established:

E(X)=p

Var(X)=c×p×(1−p)

where 0.5<c≦1.0; 0<p≦0.5. The specified global facies proportion (p)denotes the proportion of a particular facies in the facies probabilitycube, and the user-defined tuning factor or parameter controlledcontrast (c) is associated with the variance or spread of the Betadistribution. The link between the two parameters p and c and the twoparameters α and β in Beta distribution is defined as:

${\alpha = {p \times \frac{1 - c}{c}}};$$\beta = {\left( {1 - p} \right) \times \frac{1 - c}{c}}$

This equation ensures 0<α<β<1.0. As a result, the two peaks of the PDFof the Beta distribution occurs at x=0 and at x=1.0 such that the peakat x=0 is higher than the peak at x=1.0 because the target faciesproportion p is less than 0.5. If, however, the target facies proportionis more than 0.5, the complementary facies of the target facies may needto be considered to ensure that 0<α<β<1.0. Additionally, if c increasesfrom 0.5 to 1.0, the facies probability cube may tend to include morecontrast. The largest contrast is reached when c=1.0, which means thatevery value in the probability cube is either 1 or 0 and thus theprobability cube becomes a binary cube. FIG. 11 illustrates an exampleof two Probability Density Function (PDF) curves of a Beta distributionfor different parameter controlled contrasts c, which may control thespread of the Beta distribution.

In one implementation, if the global target histogram is a U-shaped betadistribution as described above, the computer application may alsoreceive the parameter controlled contrast (c) for the beta distribution.Although the global target histogram has been described as a U-shapedbeta distribution, it should be understood that in other implementationsthe global target histogram may be a different type of distribution.

At step 260, the computer application may generate a facies probabilitycube using a modified Sequential Gaussian Simulation (SGSIM) algorithm.The conventional Sequential Gaussian Simulation (SGSIM) algorithm, asdescribed in Deutsch, C. V. and Journel, A. G., 1998, GSLIB:Geostatistical Software Library and User's Guide. Oxford UniversityPress, New York, p. 369 (Deutsch and Journel, 1998), is a populargeostatistical algorithm that is used to simulate a petrophsyicalproperties distribution of any continuous variable, such as porosity orpermeability in a model (i.e., facies probability cube). It is apixel-based sequential simulation approach under a multi-Gaussianassumption.

The SGSIM algorithm first transforms well log data into standard normalvalues using a process called normal score transformation. The well logdata is transferred into standard normal values for practical purposesbecause the continuous variables may not follow a Gaussian distribution.FIG. 12 illustrates an example of a Cumulative Density Function (CDF)curve in the original probability space and a CDF curve in the normalscore space. The algorithm then proceeds by selecting a random pixel (orvoxel in 3D) in a model. Because of the multi-Gaussian assumption, theSGSIM algorithm may make a determination of a conditional cumulativedensity function (ccdf) at the selected pixel to determine theprobability of the existence of a continuous variable at the selectedpixel location. In order to make the determination of the conditionalcumulative density function (ccdf) at the selected pixel, the SGSIMalgorithm solves a local kriging system to obtain an estimated krigingmean and variance based on well data obtained from well logs that arewithin the neighborhood of the selected pixel in the model. The SGSIMalgorithm then constructs a normal distribution at the selected pixelusing the estimated kriging mean and variance. After constructing thenormal distribution at the selected pixel, the SGSIM algorithm may addthe normal distribution value to the model at the selected pixel. Thisprocess is then repeated for each pixel in the model until all of thepixels have been selected.

After applying the SGSIM algorithm to each pixel in the model in thestandard normal space, a simulation of the original (probability)variable may be obtained by performing a back transformation of thenormal score model to the original space. By using the conventionalSGSIM algorithm, the variogram model is used to control the spatialanisotropic distribution of the (probability) variable and a globaltarget histogram of the original variable is forced to be reproducedduring the back transformation.

FIG. 5 illustrates a flow diagram of a method 500 for performing amodified Sequential Gaussian Simulation (SGSIM) algorithm in accordancewith implementations of various techniques described herein. Thefollowing description of method 500 describes the process used at step260 of method 200. The modified SGSIM algorithm is a variation of theconventional SGSIM algorithm that may be used to create faciesprobability cubes such that the facies probability cubes are constrainedby the vertical facies proportion curve, the lateral proportion map, thevariogram model and the global target histogram received at steps220-250. Like the conventional SGSIM algorithm, the modified SGSIMalgorithm may first transform well log data into standard normal valuesusing a process called normal score transformation. (Step 510). Themodified SGSIM algorithm may then select a random pixel (or voxel in 3D)in a model of the facies probability cube. (Step 520). Next, in thesimulation stage, the modified SGSIM algorithm may solve a local krigingsystem at the selected pixel to determine a kriging mean m_(sk)* and akriging variance σ² _(sk) based on well data obtained from well logsthat are within the neighborhood of the selected pixel in the model.(Step 530).

Up to this step, all the procedures remain the same as the traditionalSGSIM algorithm. However, in order to constrain the facies probabilitycube to the vertical proportion curve, the modified SGSIM algorithm mayincrease or decrease the kriging mean m_(sk)* according to thedifference between a running facies proportion and a target faciesproportion at the Z-layer at selected pixel. (Step 540). The runningfacies proportion is defined as the simulated facies proportion beforethe simulation reaches the selected pixel, and the target faciesproportion is the proportion value of the target facies at the selectedpixel as indicated from the vertical proportion curve. In order tocalculate the running proportion, the computer application may perform aback transformation of the simulated pixels from the normal score valuesinto the original probability space.

In one implementation, the modified SGSIM algorithm is configured toincrease or decrease the kriging mean m_(sk)* without changing thekriging variance. As such, the chance for drawing a larger or smallerprobability value in the original space will be high unless the krigingvariance is also scaled. In order to scale the kriging mean m_(sk)*without changing the kriging variance, the modified SGSIM algorithm mayprogressively adjust the kriging mean m_(sk)* at each pixel to match thevertical proportion curve. The adjustment with regard to the verticalfacies proportion curve (VPC) is performed based on the followingformula:

$\begin{matrix}{{m_{sk}^{*{new}}\left( {i,j,k} \right)} = {{m_{sk}^{*}\left( {i,j,k} \right)} + {\frac{\lambda}{1 - \lambda}\left( \frac{{V\; P\; {C(k)}} - {p^{*}(k)}}{\sigma_{T}} \right)}}} & \left. {{Equation}\mspace{14mu} 1} \right)\end{matrix}$

where m_(sk)*(i, j, k) is the kriging mean at pixel (i, j, k), p(k)* isrunning proportion of the facies at the k-th Z-layer, VPC(k) is theproportion read at the k-th Z-layer from the VPC, λ is a servo factor in[0, 1] that indicates the correction strength, m_(sk)*^(new)(i, j, k) isthe adjusted kriging mean, and σ_(T) is the standard deviation of theglobal target histogram, i.e., σ_(T)=√{square root over (c×p×(1−p))}.

As shown in equation 1 above, if the running facies proportion at thek-th Z-layer is larger than the target vertical proportion, the krigingmean will be reduced. Alternatively, if the running facies proportion atthe k-th Z-layer is smaller than the target vertical proportion, thekriging mean will be increased. As a result, a smaller or a largernormal score value will likely be used to gear the simulated probabilityvalues in the original space towards the final facies probability cubesuch that they match the input vertical proportion curve.

If the servo factor λ is 0, the modified SGSIM algorithm does not make acorrection to the probability values. A larger servo factor (i.e.,approaching 1.0), however, corresponds with a stronger correction factorfor the reproduction of the vertical proportion curve in the resultingfacies probability cube. When the vertical proportion curve has anerratic variation cycle or trend, it may compromise the resultingprobability spatial continuity, which is typically governed by thevariogram model. Therefore, a reasonable selection of the servo factoris used to balance a trade-off between the vertical proportion curvetrend reproduction and the variogram model.

In addition to making the adjustment to the kriging mean m_(sk)* withregard to the vertical facie_(s) ^(Pr)oportion curve, the modified SGSIMalgorithm may also make an adjustment to the kriging mean m_(sk)* withregard to the lateral proportion map based on the following formula:

$\begin{matrix}{{m_{sk}^{*{new}}\left( {i,j,k} \right)} = {{m_{sk}^{*}\left( {i,j,k} \right)} + {\frac{\lambda}{1 - \lambda}\left( \frac{{V\; P\; {C(k)}} - {p^{*}(k)}}{\sigma_{T}} \right)} + {\frac{\lambda}{1 - \lambda}\left( \frac{{L\; P\; {M\left( {i,j} \right)}} - {p^{*}\left( {i,j} \right)}}{\sigma_{T}} \right)}}} & \left( {{Equation}\mspace{14mu} 2} \right)\end{matrix}$

where LPM(i, j) is lateral proportion at pixel location (i, j), whichcan be read from the LPM, and p*(i, j) is the running facies proportioncalculated from the column at (i, j). Based on equation 2, the modifiedSGSIM algorithm may constrain the facies probability cube according tothe vertical proportion curve and the lateral proportion map. (Step550).

It should be noted that two different servo factors could be used inequation 2, but by keeping the two servo factors the same the inputparameters of the algorithm may be simplified.

After scaling the kriging mean m_(sk)*, the modified SGSIM algorithm mayconstruct a normal distribution at the selected pixel using theestimated kriging mean and variance. (Step 550). After constructing thenormal distribution at the selected pixel, the modified SGSIM algorithm,like the conventional SGSIM algorithm, may add the normal distributionvalue to the model at the selected pixel. (Step 560). The modified SGSIMalgorithm may then back transform the value at the selected pixel to thetarget space and update the running proportions. (Step 570). Thisprocess is then repeated for each pixel in the model until all of thepixels have been selected.

At step 270, the computer application may perform an inverse normalscore transform in Gaussian space to the facies probability cubegenerated at step 260 in order to match the global target histogramreceived at step 250. FIG. 10 illustrates an example of a fluvial sandprobability cube generated using method 200. FIG. 15 illustrates how thefluvial sand probability cube 1530 may have been generated using the VPC1510 and the LPM 1520 as described in method 200.

In one implementation, the modified SGSIM algorithm may require thatconsistency exists between the three inputs: global facies mean (p),vertical proportion curve and lateral proportion map. The consistencymeans both the average of the vertical proportion curve and the averageof the lateral proportion map may be close to the global mean (p).Otherwise, the resulting probability cubes do not guarantee thereproduction of the input constraints.

FIG. 3 illustrates a flow diagram of a method 300 for generating amultiple facies model constrained by a facies probability cube inaccordance with implementations of various techniques described herein.The following description of method 300 is made with reference to method200 above and diagram 1500 of FIG. 15. In one implementation, method 300may be performed by a computer application.

Method 200 generates a facies probability cube by targeting one faciesat a time. As a result, the probability cubes reflect the spatialvariation of the likelihood of the selected facies. In reservoirmodeling, however, there are often more than two facies present withinthe earth. Users rarely build multiple facies probability cubes for eachfacies at one time due to the difficulty of making each faciesprobability cube consistent with each other. To overcome thisdifficulty, method 300 and method 400, described below, may be used tocreate a multiple facies model based on the facies probability cubegenerated by method 200.

At step 310, the computer application may receive a facies probabilitycube that may have been generated using method 200 described above. Thefacies probability cube as described in method 200 may also be describedas a binary (grouped) facies probability cube that includes informationpertaining to a single facies. For instance, the binary faciesprobability cube may be used to indicate the probability of whether sandexists or does not exist in an area of the earth.

At step 320, after receiving the binary facies probability cube, thecomputer application may hierarchically generate an additional binaryfacies probability cube. As such, the computer application mayrecursively generate an additional dimensions or facies on the binaryfacies cube received at step 310. In this manner, the computerapplication may repeat method 200 using the received binary faciesprobability cube. However, when repeating method 200, the computerapplication may generate a quaternary facies probability cube toindicate the presence of an additional facies. For instance, if thefacies probability cube received at step 310 represented locations wheresand exists in the survey area, at step 320, the computer applicationmay recursively evaluate the locations where sand exists in the receivedfacies probability cube such that the sand locations will be categorizedinto a different facies such as levees, crevasse or background.

Method 300 may be repeated using the newly generated binary faciesprobability cube to determine another binary facies probability cube. Assuch, method 300 may be repeated until each individual faciesprobability cube has been created. At step 325, the computer applicationmay determine whether each individual facies probability cube has beencreated. If each individual facies probability cube has been created,the computer application may proceed to step 330. If each individualfacies probability cube has not been created, the computer applicationmay return to step 310.

After creating each individual facies probability cube, at step 330, thecomputer application may receive a multi-facies training image. FIG. 13illustrates an example of a multi-facies training image.

At step 340, the computer application may generate a multiple-faciesmodel based on each individual facies probability cube. FIG. 14illustrates an example of a multiple-facies model. FIG. 15 illustrateshow the multi-facies model 1540 is generated.

FIG. 4 illustrates a flow diagram of a method 400 for generating amultiple facies model constrained by a facies probability cube inaccordance with implementations of various techniques described herein.The following description of method 400 is made with reference to method200 above and diagram 1500 of FIG. 15. In one implementation, method 400may be performed by a computer application.

At step 410, the computer application may receive a binary (grouped)facies probability cube that may have been generated using method 200described above. The facies probability cube may include the highestlevel or category of the facies group.

At step 420, the computer application may receive a multi-faciestraining image. The multi-facies training image may be used to obtain amulti-facies model using a multi-point statistics algorithm. Asmentioned above, FIG. 13 illustrates an example of a multi-faciestraining image.

At step 430, the computer application may hierarchically draw eachindividual facies at each simulated pixel in the binary faciesprobability cube. In one implementation, in order to hierarchically draweach individual facies, the computer application may use user-definedproportions for each facies of the multi-facies model. Using the definedproportions, the binary facies probability cube of step 410 and thetraining image of step 420, at step 440, the computer application maygenerate the multi-facies model using the multi-point statisticsalgorithm, described in Conditional Simulation of Complex GeologicalStructures Using Multiple Point Statistics. Mathematical Geology, v. 34,p. 1-22 (Strebelle, 2002). In one implementation, method 400 may be morepractical when there are distinct and clear facies associations in thetraining image. By generating the multi-facies model using themulti-point statistics algorithm, the computer application may generatethe multi-facies model once without using a hierarchical refiningprocess as described in method 300. As mentioned above, FIG. 15illustrates how the multi-facies model 1540 is generated.

Generally, method 300 may be used when many facies exist and one-levelfacies grouping is not enough to resolve the variability of spatialfacies heterogeneity, but it calls for many intermediate steps.Conversely, method 400 may be more practical when there are distinct andclear facies associations in the training images. As such, method 400may recursively draw facies at each simulated pixel/voxel during thesequential simulation by the following rules:

-   -   Assume there are A₁, A₂ . . . A_(M) (M>2) facies and group them        into 2 groups: A=A₁+A₂+ . . . +A_(S) and        ˜A=A_(s)+1+A_(s)+2++A_(M).    -   For a conditioning data B searched within a neighborhood of the        currently simulated pixel/voxel, read the        P(A|B)=P(A₁|B)+P(A₂|B)+ . . . +P(A₂|B) from the search tree and        combine with the underlying facies probability value (soft data)        noted as P(A|C) using the Tau model, described in Journel, A.        G., 2002, Combining Knowledge From Diverse Sources: An        Alternative to Traditional Data Independence Hypotheses.        Mathematical Geology, v. 34, p. 573-596 (Journel, 2002), to        obtain P(A|B, C).    -   Draw facies from P(A|B, C) and P(˜A|, C). If facies group A is        drawn, renormalize P(A₁), P(A₂), . . . , P(A_(S)) to 1.0 and        draw from these S facies; if facies group ˜A is drawn,        renormalize P(A_(s)+1), P(A_(s)+2), . . . , P(A_(M)) to 1.0 and        draw from these M-S facies.    -   Assign the simulated facies, which is one of A₁, A₂ . . . A_(M),        to the simulated pixel/voxel and then move to the next        simulation pixel/voxel.

FIG. 6 illustrates a computer network 600 into which implementations ofvarious technologies described herein may be implemented. The computingsystem 600 (system computer) may include one or more system computers630, which may be implemented as any conventional personal computer orserver. However, those skilled in the art will appreciate thatimplementations of various techniques described herein may be practicedin other computer system configurations, including hypertext transferprotocol (HTTP) servers, hand-held devices, multiprocessor systems,microprocessor-based or programmable consumer electronics, network PCs,minicomputers, mainframe computers, and the like.

The system computer 630 may be in communication with disk storagedevices 629, 631, and 633, which may be external hard disk storagedevices. It is contemplated that disk storage devices 629, 631, and 633are conventional hard disk drives, and as such, will be implemented byway of a local area network or by remote access. Of course, while diskstorage devices 629, 631, and 633 are illustrated as separate devices, asingle disk storage device may be used to store any and all of theprogram instructions, measurement data, and results as desired.

In one implementation, well facies logs may be stored in disk storagedevice 631. The system computer 630 may retrieve the appropriate datafrom the disk storage device 631 to predict effective permeabilitiesaccording to program instructions that correspond to implementations ofvarious techniques described herein. The program instructions may bewritten in a computer programming language, such as C++, Java and thelike. The program instructions may be stored in a computer-readablemedium, such as program disk storage device 633. Such computer-readablemedia may include computer storage media and communication media.Computer storage media may include volatile and non-volatile, andremovable and non-removable media implemented in any method ortechnology for storage of information, such as computer-readableinstructions, data structures, program modules or other data. Computerstorage media may further include RAM, ROM, erasable programmableread-only memory (EPROM), electrically erasable programmable read-onlymemory (EEPROM), flash memory or other solid state memory technology,CD-ROM, digital versatile disks (DVD), or other optical storage,magnetic cassettes, magnetic tape, magnetic disk storage or othermagnetic storage devices, or any other medium which can be used to storethe desired information and which can be accessed by the system computer630. Communication media may embody computer readable instructions, datastructures or other program modules. By way of example, and notlimitation, communication media may include wired media such as a wirednetwork or direct-wired connection, and wireless media such as acoustic,RF, infrared and other wireless media. Combinations of any of the abovemay also be included within the scope of computer readable media.

In one implementation, the system computer 630 may present outputprimarily onto graphics display 627, or alternatively via printer 628.The system computer 630 may store the results of the methods describedabove on disk storage 1029, for later use and further analysis. Thekeyboard 626 and the pointing device (e.g., a mouse, trackball, or thelike) 625 may be provided with the system computer 630 to enableinteractive operation.

The system computer 630 may be located at a data center remote from theregion were the earth formations were obtained from. The system computer630 may be in communication with the logging device described in FIG. 1(either directly or via a recording unit, not shown), to receive signalsindicating the measurements on the earth formations. These signals,after conventional formatting and other initial processing, may bestored by the system computer 630 as digital data in the disk storage631 for subsequent retrieval and processing in the manner describedabove. In one implementation, these signals and data may be sent to thesystem computer 630 directly from sensors, such as well logs and thelike. When receiving data directly from the sensors, the system computer630 may be described as part of an in-field data processing system. Inanother implementation, the system computer 630 may process seismic dataalready stored in the disk storage 631. When processing data stored inthe disk storage 631, the system computer 630 may be described as partof a remote data processing center, separate from data acquisition. Thesystem computer 630 may be configured to process data as part of thein-field data processing system, the remote data processing system or acombination thereof. While FIG. 6 illustrates the disk storage 631 asdirectly connected to the system computer 630, it is also contemplatedthat the disk storage device 631 may be accessible through a local areanetwork or by remote access. Furthermore, while disk storage devices629, 631 are illustrated as separate devices for storing input seismicdata and analysis results, the disk storage devices 629, 631 may beimplemented within a single disk drive (either together with orseparately from program disk storage device 633), or in any otherconventional manner as will be fully understood by one of skill in theart having reference to this specification.

While the foregoing is directed to implementations of varioustechnologies described herein, other and further implementations may bedevised without departing from the basic scope thereof, which may bedetermined by the claims that follow. Although the subject matter hasbeen described in language specific to structural features and/ormethodological acts, it is to be understood that the subject matterdefined in the appended claims is not necessarily limited to thespecific features or acts described above. Rather, the specific featuresand acts described above are disclosed as example forms of implementingthe claims.

What is claimed is:
 1. A method for generating one or more geologicalmodels for oil field exploration, comprising: receiving one or more wellfacies logs, a vertical facies proportion curve, a lateral proportionmap, a variogram model and a global target histogram; generating afacies probability cube using a modified Sequential Gaussian Simulation(SGSIM) algorithm, the well facies logs, the vertical facies proportioncurve, the lateral proportion map and the variogram model; matching thefacies probability cube to the global histogram; and generating thegeological models based on the matched facies probability cube.
 2. Themethod of claim 1, wherein the vertical facies proportion curve isgenerated based on an interpretation of the one or more well facieslogs.
 3. The method of claim 1, wherein the lateral proportion map isobtained by interpolating one or more facies proportions from the one ormore well facies logs.
 4. The method of claim 1, wherein the lateralproportion map is obtained by calibrating one or more seismicinterpretations of the one or more well facies logs.
 5. The method ofclaim 1, wherein the variogram model is derived based on the verticalfacies proportion curve and the lateral facies proportion map.
 6. Themethod of claim 1, wherein the variogram model is a Gaussian variogrammodel.
 7. The method of claim 1, wherein the global target histogram isa U-shaped beta distribution.
 8. The method of claim 7, whereinreceiving the global target histogram comprises receiving a parametercontrolled contrast (c) that is associated with a variance value for thebeta distribution.
 9. The method of claim 1, wherein generating thefacies probability cube using the modified Sequential GaussianSimulation (SGSIM) algorithm comprises: (a) transforming well faciesdata into a plurality of normal score values, wherein the well faciesdata is obtained from the well facies logs; (b) selecting a random pixelin a model of the facies probability cube; (c) determining a krigingmean and a kriging variance for the random pixel based on the pluralityof normal score values at the random pixel; (d) scaling the kriging meanbased on a difference between a running facies proportion, a targetfacies proportion at a Z-layer on the random pixel, the verticalproportion curve, the lateral proportion map or combinations thereof;(e) constructing a normal distribution value at the random pixel basedon the scaled kriging mean and the kriging variance; (f) adding thenormal distribution value to the random pixel in the model; and (g) backtransforming the normal distribution value at the random pixel.
 10. Themethod of claim 9, further comprising repeating steps (a)-(g) for eachpixel in the model of the facies probability cube.
 11. The method ofclaim 1, further comprising: hierarchically generating one or moreadditional facies probability cube; receiving a multi-facies trainingimage; and generating a multi-facies model based on the additionalfacies probability cubes and the multi-facies training image.
 12. Themethod of claim 1, further comprising: receiving a multi-facies trainingimage; hierarchically drawing each individual facies in the faciesprobability cube; and generating a multi-facies model using multi-pointstatistics, the multi-facies training image and each individual faciesin the facies probability cube.
 13. A computer-readable medium forgenerating one or more geological models for oil field exploration, thecomputer-readable medium having stored thereon computer-executableinstructions which, when executed by a computer, cause the computer to:receive one or more well facies logs, a vertical facies proportioncurve, a lateral proportion map, a variogram model and a global targethistogram; generate a facies probability cube using a modifiedSequential Gaussian Simulation (SGSIM) algorithm, the well facies logs,the vertical facies proportion curve, the lateral proportion map and thevariogram model; match the facies probability cube to the globalhistogram; and generate the geological models based on the matchedfacies probability cube.
 14. The computer-readable medium of claim 13,wherein the computer-executable instructions which, when executed by thecomputer, cause the computer to generate the facies probability cubeusing the modified Sequential Gaussian Simulation (SGSIM) algorithmcomprises computer-executable instructions, which, when executed by thecomputer, cause the computer to: (a) transform well facies data into aplurality of normal score values, wherein the well facies data isobtained from the well facies logs; (b) select a random pixel in a modelof the facies probability cube; (c) determine a kriging mean and akriging variance for the random pixel based on the plurality of normalscore values at the random pixel; (d) scale the kriging mean based on adifference between a running facies proportion, a target faciesproportion at a Z-layer on the random pixel, the vertical proportioncurve, the lateral proportion map or combinations thereof; (e) constructa normal distribution value at the random pixel based on the scaledkriging mean and the kriging variance; (f) add the normal distributionvalue to the random pixel in the model; and (g) back transform thenormal distribution value at the random pixel.
 15. The computer-readablemedium of claim 14, wherein the kriging mean is scaled according to:${m_{sk}^{*{new}}\left( {i,j,k} \right)} = {{m_{sk}^{*}\left( {i,j,k} \right)} + {\frac{\lambda}{1 - \lambda}\left( \frac{{V\; P\; {C(k)}} - {p^{*}(k)}}{\sigma_{T}} \right)} + {\frac{\lambda}{1 - \lambda}\left( \frac{{L\; P\; {M\left( {i,j} \right)}} - {p^{*}\left( {i,j} \right)}}{\sigma_{T}} \right)}}$where m_(sk)*(i, j, k) is the kriging mean at pixel (i, j, k), p(k)* isthe running facies proportion at the k-th Z-layer, VPC(k) is proportionread at the k-th Z-layer from the vertical proportion curve, λ is aservo factor in [0, 1] that indicates the correction strength,m_(sk)*^(new)(i, j, k) the scaled kriging mean, σ_(T) is a standarddeviation of the global histogram, i.e., σ_(T)=√{square root over(c×p×(1−p))} and LPM(i, j) is a lateral proportion at pixel location (i,j), which can be read from the lateral proportion map, and p*(i, j) is arunning facies proportion calculated from the column at (i, j).
 16. Thecomputer-readable medium of claim 14, wherein the kriging mean is scaledaccording to:${m_{sk}^{*{new}}\left( {i,j,k} \right)} = {{m_{sk}^{*}\left( {i,j,k} \right)} + {\frac{\lambda}{1 - \lambda}\left( \frac{{V\; P\; {C(k)}} - {p^{*}(k)}}{\sigma_{T}} \right)}}$where m_(sk)*(i, j, k) is a kriging mean at pixel (i, j, k), p(k)* isthe running facies proportion at the k-th Z-layer, VPC(k) is theproportion read at the k-th Z-layer from the vertical proportion curve,λ is a servo factor in [0, 1] that indicates the correction strength,m_(sk)*^(new)(i, j, k) is the scaled kriging mean, and σ_(T) is astandard deviation of the global histogram, i.e., σ_(T)=√{square rootover (c×p×(1−p))}.
 17. A system, comprising: a processor; and a memorycomprising program instructions executable by the processor to: receiveone or more well facies logs, a vertical facies proportion curve, alateral proportion map, a variogram model and a global target histogram;generate a facies probability cube using a modified Sequential GaussianSimulation (SGSIM) algorithm, the well facies logs, the vertical faciesproportion curve, the lateral proportion map and the variogram model;match the facies probability cube to the global histogram; and generateone or more geological models based on the matched facies probabilitycube, wherein the geological models are used for oil field exploration.18. The system of claim 17, wherein the global target histogram is aU-shaped beta distribution.
 19. The system of claim 17, the variogrammodel is a spherical variogram model, an exponential variogram model, aGaussian variogram model or a power variogram model.